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INTRODUCTION 


In  DME  type  instrumentation  systems,  the  ranges  from  several  distinct  locations 
to  a  common  target  are  measured  and  recorded.  If  xU^),  ytt^),  and  z(t..)  are 
the  cartesian  components  of  the  target  position  at  time  t^  in  some  reference 
coordinate  system,  and  x^,  yQ,  zq,  a=l,  H  are  the  Cartesian  components  of  the 
measurement  sites  in  the  same  reference  coordinate  system,  the  range  measure¬ 
ments  from  each  site  to  the  target  can  be  modelled  as: 

VV  =  {(x(V‘xa)2  +  (y(V_ya)2  +  (zUf)^)2}  1/2  +  eQ(ti)  a*l,M 

where  ea(t^)  is  an  error  associated  with  the  measurement  of  range  from  the  ath 
site  to  the  target.  The  process  of  estimating  tne  Cartesian  positions,  x(t^), 
y(t^),  z(t-)  from  the  measured  ranges  has  often  been  termed  multilateration. 

The  term,  trilateration,  is  a  more  common  term  used  to  denote  the  process  of 
obtaining  cartesian  positions  from  three  measured  ranges.  Several  methods  have 
been  proposed  for  multilateration,  [1  ]. 

These  methods  include  nonlinear  least  squares  and  if  M  >  4,  a  linear  least 
squares  method. 

Tne  nonlinear  least  squares  method  minimizes  the  sum  of  squares 


i  iW-  +  (y(t1- )-ya)2  +  (z(t.j)-z0)2  ): 


If  H>4,  a  linear  least  squares  method  can  be  applied  to  the  multilateration 
problem.  M-l  linear  equations  are  obtained  by  differencing  the  squares  of  the 
range  measurements.  Thus, 

Ra+l(ti)*Ro(ti)  *  2x(ti)(xa-Vl>  +  +  2z(ti>(Wl! 

+  (Wl)2  +  {Wl)2  +  {VZa+l)2 


(1) 


(2) 


(3) 


i 


> 


1 


I 


is  linear  in  the  Cartesian  positions.  Rearranging  (3)  results  in  the  M-l 
linear  equations. 

dx  x(t.)  +  d  y{tj)  +  dz  zdj)  =  P(x(ti)  i*l,  M-l 
a  y a  a 

where  d  ,  d  ,  d  are  the  direction  cosines  of  the  vector  between  measuring 
a  ya  a 

sites. 


d  «  xa”xa+l  ,  d  =  ya  yo+l  ,  d, 

a  r"“‘  y“  Za 

a.a+l  o,a+l 


=  Za-Za+1 


a,a+l 


and  R  ^  is  the  distance  between  measuring  sites 

*  <Wl>'  +  !Wl>! 

The  right  hand  side  of  (4)  is  defined  by 

Pa(ti)  *  ~  +  (*a  dx  +  ya  dy  +  dz  ) 

2Ra,a+l  ct  a  a 

The  linear  least  squares  method  estimates  x{t^ ),  y(t^),  z(t^)  by  minimizing 

V  LPa(t.)  -  (dx  ,(ti)  +  d  y(tj)  +  dz  z^))]2 

u-1  a  Ja  a 

For  applications  where  the  measuring  sites  are  all  ground  based  and  the  altitude 
of  the  vehicle  being  tracked  is  relatively  low,  both  the  linear  and  nonlinear 
least  squares  methods  experience  great  difficulty  in  obtaining  reliable  solutions. 
In  the  nonlinear  least  squares  case  this  ill-conditioning  sometimes  results  in 
the  failure  of  the  algorithm  to  converge  or  the  convergence  of  the  solution  to 
something  other  than  a  strong  local  minimum.  These  convergence  problems  exist 
even  under  the  most  ideal  conditions  of  no  error  in  the  measured  ranges.  This 
report  describes  a  nonlinear  least  squares  algorithm  for  multi lateration  which 
will  always  converge  to  a  strong  local  minimum.  In  the  ideal  case  of  no  measure¬ 
ment  error,  this  algorithm  will  converge  to  the  true  solution  provided  the 


starting  solution  is  within  certain  prescribed  limits.  The  application  of 
this  algorithm  to  the  estimation  of  a  vehicle  trajectory  from  the  measured 
ranges  of  the  General  Dynamics  RMK/MTTS  system  at  MacGregor  Range  is  described. 
The  limits  on  the  starting  solution  are  presented  for  the  case  of  perfect 
measurements  and  the  effect  of  measurement  errors  on  the  vehicle  position 
error  is  presented  for  a  wide  variety  of  geometries. 

LEAST  SQUARES  ALGORITHM  FOR  MULTI  LATE RATI ON 

The  nonlinear  least  squares  algorithm  for  multilateration  must  minimize 
M 


f(x)  •  1/2  l  (R  -f  (x))2  (2) 

a=l 

where  x  is  the  position  vector  with  coordinates  (x,  y,  z)  in  some  reference 
coordinate  system.  Rq  is  the  range  measured  by  the  ath  measurement  site  and 

fQ(x)  *  4x-XQ)2  +  (y-ya)2  +  (z-za)2  (10j 

where  (xq,  y  ,  zq)  are  the  coordinates  of  the  measuring  site  in  the  reference 
coordinate  system. 

Suppose  f(x)  is  approximated  in  the  neighborhood  of  xQ  by  a  quadratic, 

f(x)  -  f(xQ)  +  FT(xq)(x-x0)  +  1/2  (x-x0)'h(x0)(x-x0)  (11) 

where  F( o}  is  the  gradient  vector  of  f(«)»  and  H(*)  is  the  Hessian  matrix. 

Setting  the  derivative  of  f(*)  in  (11)  to  zero, 

F(x)  -  F(x0)  +  H(x0)(x-x0)  *  0  (12) 

solving  for  x  gives 

x  *  xq+H"1  (x0)F(x0)  •  (13) 

provided  H-1  (•)  exists.  The  resulting  f(x)  is 

f(x)  *  f(x0)-l/2FT  (x0)H"l(x0)F(x0)  (14) 

From  (14)  it  is  apparent  that  if  H(x0)  is  positive  definite,  f(x)<f(x0).  This 
is  the  standard  result  for  minimization  by  Newton  iteration.  If  the  Hessian 


3 


of  the  function  being  minimized  is  positive  definite,  a  sufficiently  small  step 
in  the  Newton  direction  will  always  result  in  a  decrease  of  the  objective  func¬ 


tion.  In  terms  of  the  measurement  functions  the  gradient  vector  is 


F(x)  -  -  l  F  (x) (R  -f  (x) ) , 
a*l  °  a  a 


The  Hessian  matrix  H(x)  is  given  by 

H(x)  =  l  (Fi(x)F{(x)-Ha(x)(R0-fa(x))),  (17) 

o*l 


where 

Ha(xJ  *  Vfa(x)  I3  +  (l/fa(x))3  (i-xa)(x-xa)T  (IS) 

If  a  Newton  iteration  is  used  to  solve  the  multilateration  problem,  the  Hessian 
may  often  not  be  positive  definite  so  that  the  objective  function  may  increase 
at  some  steps  and  the  direction  of  the  iteration  may  go  astray.  How  then  should 
we  proceed  at  a  point  where  the  Hessian  is  not  positive  definite?  Even  if  the 
Hessian  is  positive  definite  so  that  the  Newton  direction  is  a  descent  direction, 
the  least  squares  Newton  step  may  not  result  in  a  decrease  of  the  sum  of  squares. 
Thus,  we  have  two  problems  to  solve  in  applying  the  Newton  method  to  the  multi¬ 
lateration  problem: 

(1)  Replace  the  Hessian  matrix  with  a  matrix  H( x)  which  is  guaranteed  positive 
definite  and  a  reasonable  approximation  to  the  Hessian. 

(2)  After  obtaining  an  appropriate  H  matrix,  choose  a  step  length  o  in  the  Newton 
direction  which  results  in  a  decrease  in  the  sum  of  squares. 
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THE  GILL-MURRAY  HESSIAN  MODIFICATION 

In  order  to  insure  that  the  H  matrix  used  in  the  Newton  iteration  is  positive 
definite.  Gill  and  Murray  [2]  suggested  a  method  for  modification  of  the  Hessian 
based  on  an  LDT^  factorization  of  the  Hessian.  The  following  is  a  summary  of 
the  Gill -Murray  modification. 

If  the  Hessian  II  is  positive  definite,  then 
H  =  LDLT 

where  L  is  unit  lower  triangular  and  D  is  diagonal  with  D,.>o.  When  in  com¬ 
puting  the  L-D  factors  of  the  Hessian  we  encounter  a  D^<o,  modify  H  so  that 
H  ■  H  +  E 

so  that  H  will  be  positive  definite  and  H  =  LDLT,  D^o.  E  is  diagonal.  The 
elements  of  E  are  chosen  to  satisfy  three  requi rements ; 

(1)  All  elements  of  LD  '  are  bounded  above  by  a  chosen  value  B; 

(2)  All  elements  of  D  are  bounded  below  by  a  chosen  value  £,  and; 

(3)  H  ■  H  when  H  is  positive  definite  and  sufficiently  well  conditioned.  Thus, 

we  have  to  specify  the  parameters  B  and  £.  £  is  chosen  to  be 

S  *  max  {me| |H| | ,  me} 

where  me  is  near  the  smallest  floating  point  number  representable  by  the  com¬ 
puter  being  used  and 

I  [H||  ■  {  I 

l.j 

A  choice  of  B  which  will  satisfy  the  requirements  is 

82  *  max  {ty3,  t2,  rre}, 

where 


t2  «  max{|Hjj j } , 
j-1,3 


and 


(25) 


t i  B  max  { n j } , 
j=l,3 

with 

=  max  {|Hkj| }, 
k=j+1 ,3 


Ljp2 


D 


P 


and 


Then  define 

tj  -  ™x{|ckJ|) 
k>j 


The  elements  of  D  are  defined  as 


K  if  £>max  {|^.| ,  tj2/B2} 
kjl  if  kj|>  max  U,tj2/62} 
t j 2 / 6 2  if  t^/B2*  max  U,|i}k|} 


thus,  the  elements  of  E  are 

IE-iJ'j  it  &  max  {|^|  ,tjJ/e2} 

it  (\JijU  U,t..2/B2} 

tjVB2^-  it  max  { C » I ^ j [ > 


Thus,  we  have  the  mechanism  to  specify  a  modified  Newton  direction  which  is  a 
descent  direction.  We  must  now  specify  the  length  of  the  step  to  take  in  this 
direction. 


CHOICE  OF  STEP  LENGTH 

Having  determined  the  direction  d. 


(35) 


(27) 


(23) 


(29) 


(30) 


(31) 
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r 

i 

m  H(x0)d  =  F(x0J  (32) 

*  along  which  we  will  move  to  obtain  a  decrease  in  tnc  sum  of  squares,  we  must 

■  determine  the  lenoth  of  the  step  a  to  take  in  this  direction 

x  =  xQ  +  ad  (22) 

I  A  value  of  a*l  corresponds  to  taking  the  full  least  squares  steD  along  d.  We 

employ  the  following  strategy.  Let  the  sum  of  squares  f  along  the  direction  d 

■  be  parameterized  by  the  scalar  variable  a.  Thus,  we  denote 

n  f(a)  *  f(xQ  +  ad)  (34) 

and  the  derivative  f‘ (a)  by 

B  f’ (a)  *  F^(xq  +  ad)d  '35) 

Since  d  is  a  descent  direction,  f'(o)  =  F^(x0)d<o.  If  f(l)<f(o)  and  f'(l)<0, 

B  we  move  the  whole  least  squares  step,  ael,  so  that 

5  «  -x0  +  d  ; 

■  If  f (1 )<f (o)  and  f'(l)>o,  we  choose  a  in  the  interval  [0,1]  by  cubic  inter- 

■  polation.  In  this  case 

.  f* (1)  +  a-b 

~  0=1  -  +  2a  (37) 

*  where 

■  a  *  (b2-f‘ (l)f' (0))1/2  (38) 

b  -  3(f(0)-f(l))  +  f'(0)  +  f*(l)  (39) 

■  This  value  for  a  minimizes  f(a)  on  the  interval  [0,1]  assuming  that  f(a)  is  a 
cubic  on  this  interval.  If  f(l)>f(0),  we  halve  the  interval  considered  and 

■  repeat  the  above  strategy  including  the  cubic  interpolation  on  the  interval 

■  [0.1/2]. 

*  CONVERGENCE  CRITERIA 

■  The  Newton  iteration  is  considered  to  have  converged  when 
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(4:j 


!xi'x0  ‘  -  cl^‘xo  l+e2^»  i=  1 ’2’3 

T  i  “  _3 

for  three  consecutive  iterations.  The  constants  used  are  e-plO  and  C2=l^ 

It  is  possible  that  the  iteration  described  above  will  converge  to  a  point 
x  for  which  F(x  )=0  but  H(x  1  is  not  positive  definite.  This  situation  is 

O  O  0 

detected  when  |j  E  J | >c .  Gill  and  Murray  [2]  suggest  that  i*  this  situation 
is  encountered  the  iteration  be  continued  in  the  following  sequence.  Solve, 

LT  y  -  e,  (<1) 

where  e.  is  a  unit  vector  such  that 

J 

YEj =  nin  VEk 

k=l  ,3 

Then  choose  the  direction  for  search  as 

h  -  j  -sian(yTF(x  ))y,  if  o<||F(x  )||<c 
1  >  .if||F(«#)M-o 

APP LI  C ATI  Oil  TO  THE  RV.S/MTTS  SYSTEM 

The  current  configuration  of  the  General  Oynamics  Ri'S/t’TTS  syster  located  at 
MacGregor  Range  consists  of  eight  range  measuring  receivers.  The  geometry  of 
tiie  receiver  locations  is  given  in  the  Fic  1.  We  will  consider  the  convergence 
limits  and  position  errors  of  the  nonlinear  least  squares  algorithm  for  target 
x,y  positions  within  the  field  specified  by  the  limits  displayed  in  Fig.  1  ar  J 
for  target  z  values  of  z  =  100,  200,  300,  500,  700,  1000,  1500,  2000,  3000,  5000 
ft.  The  reference  coordinate  system  for  this  description  has  the  x,  y  plane  tan¬ 
gent  to  the  1866  Clarke  spheroid  at  C-station  with  x  is  positive  East,  y  is  pos¬ 
itive  North,  and  z  is  perpendicular  to  the  x,  y  plane  and  positive  upward.  For 
purpose  of  evaluation  we  divide  the  field  of  Fig.  1  into  cells  by  dividing  both 
the  x  and  y  limits  into  ten  equal  intervals.  Thus,  we  will  evaluate  the  conver¬ 
gence  limits  and  errors  at  121  x,  y  grid  points. 


G 


FIC.  1 


9 


jx  <\J 


We  first  consider  the  linits  on  the  starting  value  of  the  2-coordinate  of  target 
position.  Suppose  we  have  exact  measurements  of  the  range  *roir  each  o*  the  eight 


range  measuring  sites  to  the  target.  We  obtain  a  starting  solution  xpl  y0,  zo> 
for  the  nonlinear  least  squares  algorithm  by  assuring  a  zq  close  to  the  true  2 
value  and  estimating  x0,  y  ,  using  the  linear  least  souares  algorithm  described 
in  the  introduction  with  2  =  Thus,  we  are  using  the  linear  least  squares 
method  in  two  dimensions.  At  each  grid  point  on  the  field  the  nonlinear  least 


squares  algorithm  will  converge  to  a  strong  local  minimum  and  if  zQ  is  close 
enough  to  the  true  2  value  this  strong  local  minimum  will  be  the  true  value  of 
tne  target  position.  The  limits  on  the  starting  2  value  for  the  algorithm  to 
converge  to  the  true  target  position  when  presented  with  exact  range  measure¬ 
ments  are  summarized  in  the  following  table.  6zu  denotes  the  limit  on  how  far 
above  the  true  solution  that  z^  can  be  and  6z.  denotes  how  far  below  the  true 
solution  that  zQ  can  be.  The  notation  (a,  b)  in  the  table  indicates  that  the 
limit  is  somewhere  between  a  and  b  and  the  notation  >a  indicates  that  the  limit 
is  greater  than  a. 


z  (ft) 

6ZL  (ft) 

62u  (ft) 

100 

(100,150) 

(>  200) 

200 

(150,200) 

(250,300) 

300 

(200,250) 

(150,200) 

500 

(300,350) 

(150,200) 

700 

(200,250) 

(150,200) 

1000 

(0,  50) 

(200,250) 

1500 

(150,200) 

>  800 

2000 

(400,500) 

>  3000 

3000 

>  1000 

>  3000 

5000 

>  3000 

>  5000 

[«! 


O  U  O  o  o  ulu  op  Q  O 


c.  c 


Copy  available  to  DTIC  does  not 
pennit  fully  legible  reproduction 
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b  o 


C  ;C 


o  o 


X 


>763. 

0. 

0. 


*  0* 

o. 

0* 

0. 

n. 

o. 

.~  0« 

*0. 

0. 

0. 

n. 

o. 

.  0  • 

0.  " 

o. 

0. 

6« 

o. 

.  o.' 

o. 

0* 

-  3  ft  7  • 

>• 

0. ' 

«  0  • 

o. 

0. 

-535. 

-77ft. 

-73. 

-  1000 

2  •  7C0 
o 

i 

.  o. 

0. 

o. 

0. 

n. 

0. 

.  0« 

0. 

0. 

0. 

n. 

0. 

.  ‘  o. 

0. 

o. 

*0. 

o. 

o. 

.  0« 

0. 

0. 

0. 

n. 

0. 

.  “  o. 

'o. 

o. 

0. 

6* 

0.' 

.  o. 

c. 

o. 

0. 

n. 

0. 

•  C’  C-  C*  C-  C«  C'  C*  C-  O  C 


O  O  O  -3  O  O  O  • 


I 

I 

I 

I 

I 

t 

! 

i 

I 

I 

1 

I 

I 


Tiie  limits  on  tne  starting  z  solution  given  in  the  table  ere  uased  on  tiie  con¬ 
vergence  at  all  121  (x,  y)  grid  points  of  the  field.  The  limits  on  zQ  for  some 
(x,  y)  grid  points  may  be  much  greater  than  the  limits  in  the  table.  More  de¬ 
tailed  information  on  the  convergence  at  each  individual  grid  point  is  given  in 
the  preceedinr  tables.  Each  table  entry  gives  the  error  between  the  converged 
z  solution  and  the  true  2  solution.  The  heading  for  each  table  is  the  true  2 
solution  and  the  starting  solution  zQ.  The  errors  in  the  table  which  are  non- 
2ero  are  due  to  the  convergence  of  the  algorithm  to  a  local  minimum  which  is  not 
the  true  2-value. 

ESTIMATION  ERROR  DUE  TP  TRUNCATION 

After  correcting  the  measured  KTTS  ranges  for  lags  and  refraction,  one  of  the 
most  important  errors  remaining  in  the  data  is  the  two  meter  truncation  error. 
This  truncation  error  may  be  magnified  many  times  by  geometric  factors  as  will 
be  seen.  This  truncation  error  is  a  basic  limitation  of  the  MTTS  instrumenta¬ 
tion  system.  We  evaluate  the  position  error  caused  by  the  truncation  at  each  of 
the  121  (x,  y)  grid  points  previously  selected  for  each  of  the  2-coordinates, 
z  *  100,  200,  300,  500,  700,  1000,  1500,  2000,  3000,  5000  ft.  The  results  are 
presented  in  the  following  tables.  For  each  altitude  a  pair  of  error  tables  is 
given.  The  first  error  table  is  the  error  in  the  estimated  2-coordinate  for 
each  of  the  121  grid  points  in  the  field.  The  second  error  table  gives  the  error 
in  the  estimated  ground  range,  i.e.,  the  square  root  of  the  sum  of  squared  errors 
in  the  estimated  x  and  y  coordinates.  The  errors  in  both  tables  are  given  in 
feet.  The  2-coordinate  estimation  errors  are  surtmarited  in  Fic.  2  which  plots 
the  2  estimation  error  averaged  over  the  121  (x,  y)  grid  points,  the  minimum 
estimation  error  for  these  grid  points,  and  the  maximum  estimation  error  for 
these  grid  points  as  a  function  of  z.  It  is  possible  to  improve  the  estimation 
errors  due  to  the  truncation  by  making  a  simple  correction  to  the  measured  ranges. 
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If  half  of  the  truncation  error  is  added  to  each  of  the  measured  ranges,  the 
expected  value  of  the  measurement  errors  will  be  reduced  to  zero  under  reason¬ 
able  assumptions.  We  call  this  correction  a  granularity  correction.  The 
granularity  correction  for  the  RMS/MTTS  system  is  one  meter.  The  z-coordinate 
estimation  error  in  the  presence  of  a  two  meter  truncation  error  and  a  one 
meter  granularity  correction  is  summarized  In  fig.  3  where  the  average,  minimum, 
and  maximum  errors  over  the  (x,  y)  grid  points  are  plotted  as  a  ^unction  of  z. 
Comparison  of  figs.  2  and  3  show  a  significant  improvement  in  the  average  estima 
tion  error  when  the  granularity  correction  is  applied  to  the  measured  ranges. 

CONCLUSIONS  AND  FUTURE  STUDIES 

A  globally  convergent  modified  Newton  method  for  multi! ate rati  on  has  been  pre¬ 
sented.  This  algorithm  was  applied  to  the  evaluation  of  the  RMS/MTTS  measure¬ 
ment  system  at  MacGregor  Range. 

The  preceeding  tables  and  plots  show  some  severe  limitations  to  the  use  of  the 
RMS/MTTS  system  at  MacGregor  Range.  At  low  altitudes,  even  in  the  absence  of 
any  range  measurement  errors,  a  good  nominal  value  of  the  z-coordinate  of  the 
trajectory  must  be  available  in  order  that  the  trajectory  solution  converge  to 
the  true  z-solution.  The  MTTS  range  measurements  have  a  two  meter  resolution. 
The  effect  of  this  two  meter  resolution  on  the  MTTS  Cartesian  position  solution 
was  evaluated  for  a  variety  of  z-coordinates  and  ground  positions.  The  results 
of  this  evaluation  show  that  the  error  in  the  ground  plane  position  estimates 
caused  by  the  truncation  error  to  be  acceptably  small  for  all  altitudes.  How¬ 
ever,  at  low  altitudes  the  error  in  estimated  z-coordinate  is  often  unaccept¬ 
ably  large,  even  when  a  perfect  nominal  z-value  was  used  to  start  the  solution. 
It  was  demonstrated  how  a  correction  could  be  applied  to  the  range  measurements 
to  reduce  the  average  error  in  the  z-coordinate  caused  by  the  truncation  error. 


At  higher  altitudes,  say  above  3000  ft.,  the  error  in  the  estimate  of  the  z- 
component  of  position  is  of  a  more  acceptable  magnitude  and  the  modified  Newton 
algorithm  converges  to  the  true  z-solution  in  the  absence  of  measurement  error 
from  a  wide  range  of  starting  z-values. 

The  tables  presented  above  support  the  conclusion  that  the  error  in  the  esti¬ 
mated  ground  plane  position  is  acceptable  at  any  z-value.  Thus,  the  MTTS  sys¬ 
tem  offers  a  highly  desirable  measuring  system  for  estimating  the  trajectory  of 
ground  targets,  especially  if  a  terrain  map  table  look  up  is  used  to  obtain  the 
z-position  from  the  estimated  x,  y  position  estimates.  The  conclusion  that  the 
MTTS  is  attractive  for  estimating  trajectories  for  ground  targets  has  been  sup¬ 
ported  by  comparing  MTTS  estimated  position  with  positions  of  around  targets 
obtained  from  the  ALTS  laser  tracker. 

It  has  been  suggested  that  the  use  of  one  or  more  airborne  A-station  receivers 
enhances  the  estimated  z-component  of  position,  especially  at  low  altitudes. 

We  have  found  this  conclusion  to  be  correct  when  using  the  nodified  Newton 
algorithm.  The  enhancement  of  the  position  estimates  from  the  RMS /MTTS  at 
MacGregor  Range  provided  by  using  airborne  A-stations  is  the  subject  of  a  forth¬ 
coming  report. 

It  is  also  reasonable  to  suggest  that  the  position  estimate  errors  might  be 
reduced  by  using  a  position  estimation  algorithm  which  in  a  sense  averages  over 
several  x,  y  coordinates  while  estimating  a  constant  z-position.  Although  the 
computation  time  using  such  an  algorithm  would  be  significantly  greater,  this 
Increase  might  be  justified  by  the  reduced  estimation  errors.  The  use  of  such 
an  algorithm  will  be  investigated. 
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